#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use FindBin;
use lib ("$FindBin::Bin/..", "$FindBin::Bin/.", "/net/cpp-group/Leo/bin");
use chromosomes;
use linewrap;
use Data::Dump qw(dump);
use promiscuous_genes;
use db_parameters;
use ortho_species;

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_pseudo_genes


    depend on
        taxon.prot_exon_stats
        orthologs.ortholog_sets


    save
        orthologs.gene_pseudo_evidence
        orthologs.pseudo_evidence_names
		
888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

#
#	Connect to genome database
#
my $dbh = connect_to_panda();
$dbh->{RaiseError} = 0;
use get_psql_time;
my $start_time = get_psql_time($dbh);


use constant FRAME_SHIFT => 1;
use constant SINGLE_EXON => 2;
use constant DISPERSED   => 4;
use constant MULTIPLE	 => 8;
use constant PROMISCUOUS => 16;
use constant LINE1		 => 32;
use constant PSEUDO		 => 64;



#_________________________________________________________________________________________

#	get_all_genes

#_________________________________________________________________________________________
sub get_all_genes_per_species($$)
{
	my ($dbh, $species) = @_;
	print STDERR "\t\tGet a list of all genes for $species\n";
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
			DISTINCT ON (disentangled_gene_id)
            disentangled_gene_id 
		FROM 
			ortholog_sets
		WHERE 
			gene_id <> disentangled_gene_id AND
			protocol_id = $curr_protocol_id AND
			species = '$species'
        UNION
    SELECT
			DISTINCT ON (gene_id)
            gene_id
        FROM
            ens_id
		WHERE 
			species = '$species';
PL/SQLCMD
	return map{$_->[0]} @{$dbh->selectall_arrayref($sql_cmd)};
}










#_________________________________________________________________________________________

#	get_genes_with_single_exons

#		use protein loci for split genes

#_________________________________________________________________________________________
sub get_genes_with_single_exons($$\%)
{
	my ($dbh, $protocol_id, $single_exon) = @_;
	print STDERR "\t\tGet genes with single exons\n";
	for my $species($ortho_name1, $ortho_name2)
	{
		my $sql_cmd = <<"PL/SQLCMD";
		SELECT 
    			DISTINCT ON (gene_id)
    			gene_id
    		FROM 
    			taxon.gene_coding_exon_counts NATURAL JOIN 
    			ens_id
    		WHERE 
    			exon_count = 1  AND 
    			species = '$species'
        UNION
        SELECT 
                disentangled_gene_id 
            FROM 
                ortholog_sets NATURAL JOIN 
                prot_exon_stats 
            WHERE 
                real_exon_count = 1  AND 
                gene_id <> disentangled_gene_id AND
                protocol_id = $protocol_id AND
    			species = '$species'
PL/SQLCMD
        my $data = $dbh->selectall_arrayref($sql_cmd);
		print STDERR "\t\t", scalar @$data, " single exon genes for $species\n";
		++$single_exon->{$_->[0]} for (@$data);
	}
}








#_________________________________________________________________________________________

#	get_genes_with_disruptions

#		prefer proteins if orthologs

#_________________________________________________________________________________________
sub get_genes_with_disruptions($$\%\%)
{
	my ($dbh, $protocol_id, $single_disruptions, $multi_disruptions) = @_;
	print STDERR "\t\tGet genes with disruptions\n";
	print STDERR "\t\tDisruption count for all genes.\n";
	for my $species($ortho_name1, $ortho_name2)
	{
		my $sql_cmd = <<"PL/SQLCMD";
        SELECT 
                DISTINCT ON (gene_id)
                gene_id,
                MIN(disruptions_count) 
            FROM 
                taxon.prot_exon_stats NATURAL JOIN 
                ens_id 
            WHERE
                species = '$species'
            GROUP BY 
                gene_id 
            HAVING 
                MIN(disruptions_count) > 0
        UNION
        SELECT 
                disentangled_gene_id,
                disruptions_count
            FROM 
                ortholog_sets NATURAL JOIN 
                prot_exon_stats 
            WHERE 
                gene_id <> disentangled_gene_id AND
                disruptions_count > 0      AND
                protocol_id = $protocol_id AND
                species = '$species'
PL/SQLCMD
		my $cnt_multi = 0;
		my $data = $dbh->selectall_arrayref($sql_cmd) or die $dbh->errstr.".";
		for my $line (@$data)
		{
			my ($gene_id, $cnt) = @$line;
			++$single_disruptions->{$gene_id};
			if ($cnt > 1)
			{
				++$cnt_multi;
				++$multi_disruptions->{$gene_id} 
			}
		}
		print STDERR "\t\t\t", scalar @$data, " genes with disruptions from $species.\n";
		print STDERR "\t\t\t", $cnt_multi, " genes with multiple disruptions from $species.\n";
	}
}



#_________________________________________________________________________________________

#	get_dispersed_genes

#		only for orthologs

#_________________________________________________________________________________________
sub get_dispersed_genes($$\%)
{
	my ($dbh, $protocol_id, $dispersed) = @_;
	print STDERR "\t\tGet dispersed genes\n";
	for my $species($ortho_name1, $ortho_name2)
	{
		# all genes with syntenic distance are derived from ortholog_sets anyway
		my $sql_cmd = <<"PL/SQLCMD";
        SELECT 
                DISTINCT ON (disentangled_gene_id)
                disentangled_gene_id
            FROM
                orthologs.syntenic_distances NATURAL JOIN 
                orthologs.ortholog_sets
            WHERE 
                closest_neighbour_dist > 50 AND
                protocol_id = $protocol_id AND
                species = '$species';
PL/SQLCMD
		my $data = $dbh->selectall_arrayref($sql_cmd);
		++$dispersed->{$_->[0]} for (@$data);
		print STDERR "\t\t\t", scalar @$data, " dispersed from $species.\n";
	}
}



#_________________________________________________________________________________________

#	get_line1_genes

#		only for orthologs

#_________________________________________________________________________________________
sub get_line1_genes($$\%)
{
	my ($dbh, $protocol_id, $line1) = @_;
	print STDERR "\t\tGet L1 transposable elements genes\n";
	for my $species($ortho_name1, $ortho_name2)
	{
		# select all gene ids from any ortholog set where this is an interpro line1 hit
		# for either species
		my $sql_cmd = <<"PL/SQLCMD";
        SELECT 
                DISTINCT ON (disentangled_gene_id)
                disentangled_gene_id
            FROM 
                orthologs.ortholog_sets NATURAL JOIN
                (
                    SELECT 
                           ortholog_id,
                           protocol_id
                        FROM 
                            prot_interpro_matches NATURAL JOIN 
                            ens_id JOIN 
                            orthologs.ortholog_sets USING (gene_id, species)
                        WHERE 
                            interpro_ac = 'IPR004244' AND
                            protocol_id = $protocol_id
                ) AS a
            WHERE 
                species = '$species';
PL/SQLCMD
		my $data = $dbh->selectall_arrayref($sql_cmd);
		++$line1->{$_->[0]} for (@$data);
		print STDERR "\t\t\t", scalar @$data, " L1 transposable elements from $species.\n";
	}
	
}

#sub get_promiscuous_genes($$$)


#_________________________________________________________________________________________

#	get_ortholog_sets

#_________________________________________________________________________________________
sub get_ortholog_sets($$\%\%\%)
{
	my ($dbh, $protocol_id, $ortholog_sets, $ortholog_counts, $gene_id_to_ortholog_id) = @_;
	print STDERR "\t\tGet sets of ortholog predictions\n";
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
		disentangled_gene_id,
        ortholog_id
    FROM
        orthologs.ortholog_sets
    WHERE 
        protocol_id = $protocol_id;
PL/SQLCMD
	# sort into species
	for my $line (@{$dbh->selectall_arrayref($sql_cmd)})
	{
		my ($gene_id, $ortholog_id) = @$line;
		my $species = ($gene_id =~ /$ortho_gene1/) ? 0 : 1;
		$ortholog_sets->{$ortholog_id}{$species}{$gene_id}=undef;
		$gene_id_to_ortholog_id->{$gene_id} = $ortholog_id;
	}

	# count how many each
	for my $ortholog_id(keys %$ortholog_sets)
	{
		$ortholog_counts->{$ortholog_id}[0] = 
										scalar keys %{$ortholog_sets->{$ortholog_id}{0}};
		$ortholog_counts->{$ortholog_id}[1] = 
										scalar keys %{$ortholog_sets->{$ortholog_id}{1}};
	}
}







#-----------------------------------------------------------------------------------------
#

#	save names for pseudos

#		names for the different types of evidence for retro pseudo genes
#
#-----------------------------------------------------------------------------------------
sub save_pseudo_evidence_names($)
{
	#
	#	set up the names for the different reasons for pseudo genes
	#
	my %pseudo_types = (
					  "Single-exon"                 					=> 2,
					  "Dispersed"                   					=> 4,
					  "Multiple Frame-shifts"							=> 8,
					  "Promiscuous"										=> 16,
					  "L1 transposable element"							=> 32,
	);
	my %pseudo_ids =reverse %pseudo_types;
	my @reasons_names;
	for (my $i = 0; $i <= 63; ++$i)
	{
		my @type_names;
		for my $pseudo_id(keys %pseudo_ids)
		{
			if ($i & $pseudo_id)
			{
				push(@type_names, $pseudo_ids{$pseudo_id});
			}
		}
		if (!($i & 8) && ($i & 1))
		{
			push(@type_names, "Frame-shift");
		}
		my $reasons_name =join (" and ",  @type_names);
		if ($i >= 5)
		{
			$reasons_names[$i] =
			$reasons_names[$i + 64] = "Pseudo: ".$reasons_name;
		}
		else
		{
			$reasons_names[$i] = $reasons_name;
			$reasons_names[$i + 64] = "Error: This should not be labelled as a pseudogene";
		}
	}
	
	
	

	print STDERR "\tSave names for component evidence for pseudos\n";
	$dbh->do("TRUNCATE orthologs.pseudo_evidence_names");
	$dbh->do("COPY orthologs.pseudo_evidence_names".
			"(reason, reason_id) FROM STDIN");
	$dbh->pg_putline("$reasons_names[$_]\t$_\n") for (0 .. $#reasons_names);
	$dbh->pg_endcopy();
	$dbh->do("ANALYSE orthologs.pseudo_evidence_names;");
}





#-----------------------------------------------------------------------------------------
#
#	integrate_pseudo_evidence_per_gene

#		For genes, check for the component evidence for pseudo-ness, i.e.
#			dispersal, disruptions, single exon

#
#-----------------------------------------------------------------------------------------
sub integrate_pseudo_evidence_per_gene(\@\%\%\%\%\%\%\%\%)
{
	my ($genes, $single_exon, $single_disruptions, 
				$multi_disruptions, $dispersed, 
				$line1, $promiscuous_genes,
				$pseudo_evidence_per_gene, $pseudo_genes) = @_;
	print STDERR ("\t\tDiagnosing pseudos...\n");
	# for each gene
	for my $gene_id(@$genes)
	{
		my $pseudo_evidence = 0;

		#
		# pseudo 
		#
		
			# if multiple disruptions
			if (exists $multi_disruptions->{$gene_id})
			{
				$pseudo_evidence |= MULTIPLE + 0 + PSEUDO;
			}
			# line 1
			if (exists $line1->{$gene_id})
			{
				$pseudo_evidence |= LINE1 + 0 + PSEUDO;
			}

			# promiscuous gene
			if (exists $promiscuous_genes->{$gene_id})
			{
				$pseudo_evidence |= PROMISCUOUS  + 0 + PSEUDO;
			}

		#
		# contributes toward pseudo 
		#
			# single exon
			if (exists $single_exon->{$gene_id})
			{
				$pseudo_evidence |= SINGLE_EXON;
			}
			# *any* disruption
			if (exists $single_disruptions->{$gene_id})
			{
				$pseudo_evidence |= FRAME_SHIFT;
			}
			
			# dispersed
			if (exists $dispersed->{$gene_id})
			{
				# pseudo if dispersed and single exon or frameshift
				$pseudo_evidence |= DISPERSED;
				if ($pseudo_evidence & (SINGLE_EXON | FRAME_SHIFT) )
				{
					$pseudo_evidence |= PSEUDO;
				}
			}

		# save
		if (($pseudo_evidence & PSEUDO) != 0)
		{
			$pseudo_genes->{$gene_id}++;
		}
		$pseudo_evidence_per_gene->{$gene_id} = $pseudo_evidence;
	}
}




#_________________________________________________________________________________________

#	set_pseudo_types_for_orthologs

#_________________________________________________________________________________________
sub set_pseudo_types_for_orthologs(\%\%\%\%)
{
	my ($orthologs, $pseudos, $ortholog_counts, $gene_id_to_ortho_type) = @_;
	print STDERR ("\t\tAssigning orthologs to pseudos or orphans...\n");
	for my $ortholog_id (keys %$orthologs)
	{
		my @pseudo_counts = (0,0);
		for my $species(0,1)
		{
			my $ortholog_set = $orthologs->{$ortholog_id}{$species};
			for my $gene_id(keys %$ortholog_set)
			{
				if (exists $pseudos->{$gene_id})
				{
					$gene_id_to_ortho_type->{$gene_id} = 'pseudo';
					++$pseudo_counts[$species];
				}
			}
		}

		#
		#	Set orphaned by pseudo
		#
		
		#		ignore Orphan inparalogs
        next unless (	$ortholog_counts->{$ortholog_id}[0] &&
						$ortholog_counts->{$ortholog_id}[1]);

		# Orphaned by pseudo if all orthologs in one species are pseudos
		next unless ($ortholog_counts->{$ortholog_id}[0] == $pseudo_counts[0] ||
					 $ortholog_counts->{$ortholog_id}[1] == $pseudo_counts[1]);
		
		for my $species(0,1)
		{
			my $ortholog_set = $orthologs->{$ortholog_id}{$species};
			for my $gene_id(keys %$ortholog_set)
			{
				# orphaned unless already pseudo
				$gene_id_to_ortho_type->{$gene_id} = 'Orphan by pseudo'
					unless (exists $gene_id_to_ortho_type->{$gene_id});
			}
		}
	}
}


#_________________________________________________________________________________________

#	get_counts_for_pseudo_evidence_type_for_ones_and_manys

#_________________________________________________________________________________________
sub get_counts_for_pseudo_evidence_type_for_ones_and_manys(\%\%\%\%)
{
	my ($orthologs, $ortholog_counts, $pseudo_evidence_per_gene, $count_pseudo_evidence_genes) = @_;
	print STDERR ("\t\tGet counts for genes matching the different pseudo evidence categories...\n");
	for my $ortholog_id (keys %$orthologs)
	{
		my $ortholog_type = "many";
		if ($ortholog_counts->{$ortholog_id}[0] == 0 ||
			$ortholog_counts->{$ortholog_id}[1] == 0)
		{
			$ortholog_type = "Orphan paralog";
		}
		elsif (	$ortholog_counts->{$ortholog_id}[0] +
				$ortholog_counts->{$ortholog_id}[1] == 2)
		{		
			$ortholog_type = "1 to 1"
		}
		for my $species(0,1)
		{
			my $ortholog_set = $orthologs->{$ortholog_id}{$species};
			for my $gene_id(keys %$ortholog_set)
			{
				# count the number of genes which belong to either 1 to 1 or many
				#	and which have this type of pseudo evidence
				++$count_pseudo_evidence_genes->{$pseudo_evidence_per_gene->{$gene_id}}
							{$ortholog_type}[$species];
			}
		}
	}
	#print STDERR dump($count_pseudo_evidence_genes);
}




#_________________________________________________________________________________________

#	get_pseudo_summary

#_________________________________________________________________________________________
sub get_pseudo_summary(\%\%)
{
	my ($pseudo_evidence_gene_counts, $pseudo_summary) = @_;
	print STDERR ("\t\tSummarise pseudo/non-pseudo counts...\n");

	# types of things we want to save:
	# Things which are supersets of the below
	my @pseudo_evidence_types	= (	FRAME_SHIFT,
									SINGLE_EXON,
									DISPERSED,
									DISPERSED + 0 + FRAME_SHIFT,
									DISPERSED + 0 + SINGLE_EXON,
									MULTIPLE,
									LINE1,
									PROMISCUOUS,
									PSEUDO);
	# Evidence which exactly match the below
	my @non_pseudo_variants		= (DISPERSED, FRAME_SHIFT, SINGLE_EXON);
	my %pseudo_variants_names;
	@pseudo_variants_names{@pseudo_evidence_types} = 
	(   
		"10 All Frame-shifts",
		"20 All Single-exons",
		"30 All Dispersed",
		"40 Pseudos: Dispersed and Frame-shift",
		"50 Pseudos: Dispersed and Single-exon",
		"60 Pseudos: Multiple Frame-shifts",
		"70 Pseudos: L1 transposable elements",
		"80 Pseudos: Promiscuous gene family fragments",
		"90 Pseudos");
	$pseudo_variants_names{DISPERSED + 100}   	= "91 Non-pseudo Dispersed";
	$pseudo_variants_names{FRAME_SHIFT + 100} 	= "92 Non-pseudo Frameshifts";
	$pseudo_variants_names{SINGLE_EXON + 100}	= "93 Non-pseudo Single Exons";

	#
	#	start with zero counts
	#
	for my $many_or_1("1 to 1", "many", "Orphan paralog")
	{
		for my $species(0,1)
		{
			for my $type_name(values %pseudo_variants_names)
			{
				$pseudo_summary->{$type_name}{$many_or_1}[$species] = 0;
			}
		}
	}

	#
	#	Fit gathered pseudo evidence into the above categories
	#
	for my $reason(keys %$pseudo_evidence_gene_counts)
	{
		for my $many_or_1("1 to 1", "many", "Orphan paralog")
		{
			for my $species(0,1)
			{
				next unless defined  $pseudo_evidence_gene_counts->{$reason}{$many_or_1}[$species];


				# if type is superset, add count to type
				for my $type(@pseudo_evidence_types)
				{
					next unless (($reason & $type) == $type);
					my $type_name = $pseudo_variants_names{$type};
					$pseudo_summary->{$type_name}{$many_or_1}[$species] += 
							$pseudo_evidence_gene_counts->{$reason}{$many_or_1}[$species];

							
				}
	
				# if is type but non-pseudo, note that
				for my $type(@non_pseudo_variants)
				{
					next if ($reason & $type) == 0;
					next unless ($reason & PSEUDO) == 0 ;
					my $type_name = $pseudo_variants_names{$type + 100};
					$pseudo_summary->{$type_name}{$many_or_1}[$species] += 
							$pseudo_evidence_gene_counts->{$reason}{$many_or_1}[$species];
				}
			}
		}
	}

}




#_________________________________________________________________________________________

#	save_pseudo_summary_to_db

#_________________________________________________________________________________________
sub save_pseudo_summary_to_db($$\%)
{
	my ($dbh, $protocol_id, $pseudo_summary) = @_;
	print STDERR ("\t\tSave pseudo/non-pseudo counts to panda...\n");

	$dbh->do("DELETE FROM orthologs.pseudo_analysis_summary where protocol_id = $protocol_id;");
	$dbh->do("VACUUM ANALYSE orthologs.pseudo_analysis_summary;");
	$dbh->do("COPY orthologs.pseudo_analysis_summary ".
			 "(protocol_id, category, species, ortholog_type, gene_count) FROM STDIN") or die $dbh->errstr.".";

	my @species_names = ($ortho_name1, $ortho_name2);
	for my $summary_key(sort keys %$pseudo_summary)
	{
		my $category = $summary_key;
#		$category =~ tr/0-9//d;
		for my $species(0, 1)
		{
			for my $many_or_1("1 to 1", "many", "Orphan paralog")
			{
				$dbh->pg_putline(join("\t",  
									  $protocol_id,
									  $category,
									  $species_names[$species],
									  $many_or_1, 
									  $pseudo_summary->{$summary_key}{$many_or_1}[$species]."\n"));
			}
		}
	}
	$dbh->pg_endcopy();
}











#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#			Main logic starts here


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888





	#
	#	Each pseudo evidence, e.g. single exon, has an ID and name
	#
	print STDERR ("\tSetting retro names in panda...\n");
	save_pseudo_evidence_names($dbh);
	
	#
	#	List of species
	#
	my @ortho_species = ($ortho_name1, $ortho_name2);
	
	#
	#	Prepare panda tables
	#
	$dbh->do("DELETE FROM orthologs.gene_pseudo_evidence where protocol_id = $curr_protocol_id;");
	$dbh->do("VACUUM ANALYSE orthologs.gene_pseudo_evidence;");


	# get protocol_name
	my $protocol_name = get_ortholog_name($curr_protocol_id);
	print STDERR "\tOrthology protocol [$protocol_name]...\n";


	
	#
	#	Get all the genes
	#
	my (@genes1,  @genes2,  @all_genes);
	push(@genes1, get_all_genes_per_species($dbh, $ortho_name1));
	push(@genes2, get_all_genes_per_species($dbh, $ortho_name2));
	push(@all_genes, @genes1,  @genes2);
	

	#
	#	Get components evidence for pseudos
	#
	my (%single_exon,	%single_disruptions, %multi_disruptions, %dispersed);
	my (%line1, %promiscuous_genes);
	get_genes_with_single_exons($dbh, $curr_protocol_id, 
								%single_exon);
	get_genes_with_disruptions($dbh, $curr_protocol_id, 
								%single_disruptions, %multi_disruptions);
	get_dispersed_genes($dbh, $curr_protocol_id, %dispersed);
	get_line1_genes($dbh, $curr_protocol_id, %line1);
	get_promiscuous_genes($dbh, $curr_protocol_id, %promiscuous_genes, $ortho_name1, $ortho_name2);



	#
	#	Get data from orthologs.ortholog_sets
	#
	my (%ortholog_sets, %ortholog_counts, %gene_id_to_ortholog_id);
	get_ortholog_sets($dbh, $curr_protocol_id, %ortholog_sets, %ortholog_counts, %gene_id_to_ortholog_id);




	#
	#	Calculate pseudos integrating evidence per gene
	#
	my (%pseudo_genes, %pseudo_evidence_per_gene);
	integrate_pseudo_evidence_per_gene(	@all_genes, 
										# component evidence
										%single_exon, %single_disruptions, 
										%multi_disruptions, %dispersed, 
										%line1, %promiscuous_genes,
										
										# integrated evidence
										%pseudo_evidence_per_gene, 
										
										# retrieve pseudo and real genes
										%pseudo_genes);





	#
	#	count how many genes are dispersed, have disruptions etc.
	#		Use "naive" ortholog counts before pseudos are taken out
	#
	my (%count_pseudo_evidence_genes);
	get_counts_for_pseudo_evidence_type_for_ones_and_manys(%ortholog_sets, 
														   %ortholog_counts, 
														   %pseudo_evidence_per_gene, 
														   %count_pseudo_evidence_genes);


	#
	#	sort into categories for printing out
	#
	my %pseudo_summary;
	get_pseudo_summary(%count_pseudo_evidence_genes, %pseudo_summary);
	save_pseudo_summary_to_db($dbh, $curr_protocol_id, %pseudo_summary);



	#
	#	mark pseudos among orthologs:
	#
	my (%gene_id_to_ortho_type);
	set_pseudo_types_for_orthologs(	%ortholog_sets, %pseudo_genes, %ortholog_counts,
								    %gene_id_to_ortho_type);




	#
	#	Enter into panda ortholog_sets and pseudo evidence
	#
	{
		print STDERR "\t\tMark pseudos among orthologs.ortholog_sets...\n";
		$dbh->do("CREATE TEMP TABLE t_ortholog_sets(disentangled_gene_id TEXT, ortholog_type TEXT);") or die $dbh->errstr.".";
		$dbh->do("COPY t_ortholog_sets( ".
				 "disentangled_gene_id, ortholog_type) FROM STDIN") or die $dbh->errstr.".";
		for my $gene_id(keys %gene_id_to_ortho_type)
		{
			die unless exists $gene_id_to_ortholog_id{$gene_id};
			$dbh->pg_putline("$gene_id\t$gene_id_to_ortho_type{$gene_id}\n");
		}
		$dbh->pg_endcopy();
	
		my $sql_cmd = <<"PL/SQLCMD";
		UPDATE
				ortholog_sets
			SET 
				ortholog_type = t_ortholog_sets.ortholog_type
			FROM 
				t_ortholog_sets
			WHERE
				protocol_id = $curr_protocol_id AND
				ortholog_sets.disentangled_gene_id = t_ortholog_sets.disentangled_gene_id;
		VACUUM ANALYSE orthologs.ortholog_sets;
PL/SQLCMD
		$dbh->do($sql_cmd) or die $dbh->errstr.".";
	}


	#
	#	Component evidence for pseudo genes
	#
	{
		print STDERR "\t\tEnter Component evidence for pseudo genes into panda...\n";
		$dbh->do("COPY orthologs.gene_pseudo_evidence(".
						"disentangled_gene_id, gene_id, reason_id, protocol_id ) FROM STDIN") or die $dbh->errstr.".";
		$dbh->pg_putline("$_\t$_\t$pseudo_evidence_per_gene{$_}\t$curr_protocol_id\n") for (@all_genes);
		$dbh->pg_endcopy();
	
		print STDERR "\t\tAnalysing tables panda...\n";
		$dbh->do("ANALYSE orthologs.gene_pseudo_evidence;");
	}
    


#_________________________________________________________________________________________
#	
#		orthologs.gene_pseudo_stats
#
#
    print STDERR "\t\tUpdating gene_pseudo_stats...\n";
    $dbh->do( <<"PL/SQLCMD");
    DELETE FROM taxon.gene_pseudo_stats WHERE protocol_id = $curr_protocol_id;
    INSERT INTO taxon.gene_pseudo_stats
        (gene_id, max_exon_count, min_disruptions_count,protocol_id) 
        SELECT 
                DISTINCT ON (gene_id) 
                gene_id, 
                MAX(real_exon_count) AS max_exon_count, 
                MIN(disruptions_count) AS min_disruptions_count,
                $curr_protocol_id
            FROM 
                taxon.prot_exon_stats RIGHT JOIN 
                ens_id USING (prot_id)
                GROUP BY gene_id;
    UPDATE taxon.gene_pseudo_stats
        SET dispersed = syntenic_distances.closest_neighbour_dist < 50
        FROM orthologs.syntenic_distances
        WHERE
            syntenic_distances.gene_id = gene_pseudo_stats.gene_id AND
            syntenic_distances.protocol_id = gene_pseudo_stats.protocol_id;
    VACUUM ANALYSE taxon.gene_pseudo_stats;
PL/SQLCMD


#_________________________________________________________________________________________
#	
#		taxon.gene_pseudo_stats
#			redundant here if pipeline run sequentially
#
#
print STDERR "\tUpdate taxon.gene_pseudo_stats\n";
$dbh->do( <<"PL/SQLCMD") or die $dbh->errstr.".";
DELETE FROM taxon.gene_pseudo_stats WHERE protocol_id = $curr_protocol_id;
INSERT INTO taxon.gene_pseudo_stats
    (gene_id, max_exon_count, min_disruptions_count,protocol_id) 
    SELECT 
            DISTINCT ON (gene_id) 
            gene_id, 
            MAX(real_exon_count) AS max_exon_count, 
            MIN(disruptions_count) AS min_disruptions_count,
            $curr_protocol_id
        FROM 
            taxon.prot_exon_stats RIGHT JOIN 
            ens_id USING (prot_id)
            GROUP BY gene_id;
UPDATE taxon.gene_pseudo_stats
    SET dispersed = syntenic_distances.closest_neighbour_dist < 50
    FROM orthologs.syntenic_distances
    WHERE
        syntenic_distances.gene_id = gene_pseudo_stats.gene_id AND
        syntenic_distances.protocol_id = gene_pseudo_stats.protocol_id;
VACUUM ANALYSE taxon.gene_pseudo_stats;
PL/SQLCMD



log_pipeline_stage($dbh, $curr_protocol_id, 37, 'table pseudo_genes', $start_time);
$dbh->disconnect();



 
